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Abstract 

This paper presents a micromechanical study of unsaturated granular media in 
the pendular regime, based upon numerical experiments using the discrete element 
method, compared to a microstructural elastoplastic model. Water effects are taken 
into account by adding capillary menisci at contacts and their consequences in terms 
of force and water volume are studied. Simulations of triaxial compression tests are 
used to investigate both macro and micro-effects of a partial saturation. The results 
provided by the two methods appear to be in good agreement, reproducing the 
major trends of a partially saturated granular assembly, such as the increase in 
the shear strength and the hardening with suction. Moreover, a capillary stress 
tensor is exhibited from capillary forces by using homogenisation techniques. Both 
macroscopic and microscopic considerations emphasize an induced anisotropy of 
the capillary stress tensor in relation with the pore fluid distribution inside the 
material. In so far as the tensorial nature of this fluid fabric implies shear effects 
on the solid phase associated with suction, a comparison has been made with the 
standard equivalent pore pressure assumption. It is shown that water effects induce 
microstrural phenomena that cannot be considered at the macro level, particularly 
when dealing with material history. Thus, the study points out that unsaturated soil 
stress definitions should include, besides the macroscopic stresses such as the total 
stress, the microscopic interparticle stresses such as the ones resulting from capillary 
forces, in order to interpret more precisely the implications of the pore fluid on the 
mechanical behaviour of granular materials, [published, DOI: 10.1002/nag.767] 

Key words: micromechanics; granular materials; unsaturated; DEM; capillary 
forces; microstructure 
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1 INTRODUCTION 



Macroscopic properties of granular materials such as soils depend 
on particle interactions. In dry granular materials, interparticle 
forces are related to the applied external loads as different stud- 
ies have shown In unsaturated soils subjected to capillary 
effects, new features must be accounted for in order to under- 
stand properly their behaviour. The presence of water leads to 
the formation of water menisci between neighboring grains, intro- 
ducing new interparticle forces. The effects of these forces depend 
on the degree of saturation of the medium. For low water con- 
tent level corresponding to disconnected liquid bridges between 
grains, capillary theory allows the force induced by those bridges 
to be linked to the local geometry of the grains and to the matric 
suction or capillary pressure inside the medium [3] . Since the dis- 
connected menisci assumption is not valid for high water content 
levels due to water percolation, we consider here only the unsat- 
urated state where the discontinuity of the water phase can be 
assumed, the so-called pendular regime. 

There has been a wide debate on the various interpretations for 
the mechanical behaviour of unsaturated soils. At early stages of 
soil mechanics, Terzaghi [1] first introduced the concept of effec- 
tive stress for the particular case of saturated soils enabling the 
conversion of a multiphase porous medium into a mechanically 
equivalent single-phase continuum. In unsaturated soils, water in- 
duced stresses are still debated. The common practice PjB] is to 
use the suction or a modified version as a second stress variable 
within a complete hydro-mechanical framework. An alternative 
method is introduced to develop homogenisation techniques in or- 
der to derive stress-strain relationships from forces and displace- 
ments at the particle level as proposed in |[T] for dry granular 
materials. The basic idea is to consider the material as repre- 
sented by a set of micro-systems, postulating that the behaviour 
of a material volume element depends on the intergranular inter- 
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actions belonging to this volume. We propose here to extend this 
micro-mechanical approach to unsaturated granular materials as 
proposed by Li [7], Jiang et al. |B] or Lu and Likos p. 

Along these lines we present two micromechanical models which 
take into account capillary forces. The first one is a three dimen- 
sional numerical model based on the Discrete Element Method 
(hereafter designed as the DEM model) pioneered by Cundall 
and Strack ||2], and the second one is an analytical model (here- 
after designed as the microstructural model) recently proposed 
by Hicher and Chang fW\ . 

The microstructural model is a stress-strain relation which con- 
siders interparticle forces and displacements. Thanks to analyt- 
ical homogenisation/localisation techniques, the macroscopic re- 
lation between stress and strain can be derived. 

In the DEM model, a granular medium is modelled by a set of 
grains interacting according to elementary laws. Direct simula- 
tions are carried out on grain assemblies, computing the response 
of the material along a given loading path. 

By studying their effects under triaxial loading, we investigate 
capillary forces implications at the macroscopic level, and offer an 
insight into the unsaturated soil stress framework by introducing 
a capillary stress tensor as a result of homogenization techniques. 

2 UNSATURATED SOIL STRESSES 

2.1 Macroscopic views 

Macroscopic interpretations of the mechanical behaviour of un- 
saturated soils have been mainly developed in the framework of 
elastoplasticity |ll] . Most of these models consider that the strain 
tensor is governed by the net stress tensor cr^j — (wa being 
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the pore air pressure) and the matric suction or capillary pressure 
Ua — Uw {uw being the pore water pressure) inside the medium 
[T2lfT3] . In particular, they consider a new yield surface, called 
Loading Collapse (LC) surface in the plane {{(Jij—Ua5ij) ^{ua—Uw)) 
which controls the volume changes due to the evolution of the de- 
gree of saturation for a given loading path. As a matter of fact, all 
these formulations can be considered as extensions of the relation- 
ship initially proposed by Bishop and Blight [[H] for unsaturated 
soils: 

where xi^r) is called the effective stress parameter or Bishop's 
parameter, and is a function of the degree of saturation Sr of the 
medium (x = for a dry material, x = 1 for a fully saturated 
material). 

Obviously, since the effective stress principle is by definition a 
macroscopic concept, several authors (Lu and Likos [S] or Li [7]) 
have proposed to use a micromechanical approach for the effec- 
tive stress principle. In order to further study this micromechan- 
ical approach to study unsaturated soil stresses, we propose here 
a micromechanical analysis of the problem, examining the local 
water induced effects through a set of simulated laboratory ex- 
periments. 

2.2 Micromechanical interpretation 

Let us consider a Representative Volume Element (RVE) of a 
wet granular material, subjected to an assigned external loading. 
When the water content decreases inside a saturated granular 
sample, the air breaks through at a given state. The capillary 
pressure {ua — u^) corresponding to that point is called the air- 
entry pressure, and strongly depends on the pore sizes. There- 
after, the sample becomes unsaturated and capillary forces start 
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to grow due to interface actions between air and water. Since 
the the gaseous phase is discontinue, this is the capillary regime. 
From this state, a constant decrease in the degree of saturation 
corresponds to a gentle increase in pore water pressure. The pen- 
dular regime starts when the water phase is no longer continuous. 
In this state, fluids equilibrium is obtained by the vapor pressure. 
Analytical and experimental results demonstrate that cap- 

illary effects at particle contacts produces a kind of bond between 
particles as a result of menisci (FigH]). Liquid bridges may form 
between some pairs of adjoining particles not necessarily in con- 
tact, generating an attractive capillary force between the bonded 
particles. If the drying process continues, these water bridges be- 
gin to fail, starting from the non-contacting grains, until the com- 
plete disappearance of capillary forces inside the assembly. 

As the pendular regime is considered throughout this paper, wa- 
ter is considered to be solely composed of capillary menisci: each 
liquid bridge is assumed to connect only two particles. There- 
fore, two types of forces coexist within the granular medium. 
For dry contacts, a contact force develops between contacting 
granules. This repulsive force, that is a function of the relative 
motion between the contacting grains, is usually well described 
by an elastoplastic contact model. For water bonded particles, a 
specific attractive force exists. This water-induced attractive in- 
teraction can be described by a resulting capillary force, rather 
than by a stress distribution as mentionned by Haines |[T5] or 
Fisher |j3]. This capillary force is a function of the bridge volume, 
of the size of particles, and of the fluid nature (see section 3.1.1 
for the details). The objective of this section is to derive, in a 
simple manner, an expression relating the overall stress tensor 
within the RVE to this internal force distribution. 

For this purpose, the Love static homogenisation relation 
is used. This relation expresses the mean stress tensor a within 
a granular volume ^ as a function of the external forces F^^^'P 
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applied to the particles p belonging to the boundary dV of the 
volume: 

= r,'^ Fr'-'x-; (2) 

^ pedV 



where Xj are the coordinates of the particle p with respect to a 
suitable frame. It is worth noting that this relation is valid what- 
ever the nature of the interactions between grains. 



Taking into account the mechanical balance of each particle of the 
volume V (including the boundary dV), Eq.(2) can be written 
as: 

1 N N 

<^., = T7 s E prtf (3) 

^ V=lq=l 



where N is the number of particles within the volume, F^'^ is 
the interaction force exerted by the particle q onto the particle j9, 
and /^'^ is the branch vector pointing from particle q to particle 



As we consider partially saturated granular media, two indepen- 
dent kinds of interparticle forces can be distinguished: 

(i) if particles p and q are in contact, a contact force F^^^ exists. 

(ii) if particles p and q are bonded by a liquid bridge, a capillary 
force F^^p exists. 

Actually, depending on the local geometry, a liquid bond can ex- 
ist between two grains in contact. In that case, solid contacts 
are surrounded by the continuous liquid phase providing the si- 
multaneity of contact and capillary forces. The two contributions 
have therefore to be accounted for by summation. 
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Finally, in all cases and for any couple {p, q)e[l, N]'^, it can be 
written that: 

F'^^' = FZt + FX (4) 
Thus, by combining Eqs.(2) and (4), it follows that: 

N N ■] N N 

^ij ~ T/ ^contyj I -rr ^cap,i^j 

^ p=lq=l ^ 15=1 <?=! 



As a consequence, Eq.(5) indicates that the stress tensor is split 
into two components: 

_ cont I cap /n\ 

= (^^J + (6) 

- A first component afj^^ = y E^i E^i ^cont,^^!'^ accounting 
for the contact forces transmitted along the contact network. 

- A second component a'^j^ = y E^i E^i F^^^lf^ represent- 
ing the capillary forces existing within the assembly. 

It is to be noted that crfj"* is a stress quantity standing for in- 
tergranular contact forces in the same way as in saturated or 
dry conditions. Considering the concept as initially introduced 
by Terzaghi, crfj"* plays the role of the so-called effective stress 
by governing soil deformation and failure. Besides, a^j^ is the ten- 
sorial attribute to capillary water effects or suction, by extension. 
By analogy with Eq.(l), we can therefore define a microstructural 
effective stress 

^cont ^ cap (ry\ 



where cr^j could be affiliated to net stress, representing the ap- 
parent stress state in the material. Compared with Eq.(l), where 
the effect of water is intrinsically isotropic, a'^j^ implies a tenso- 
rial attribute to the water effects. 



In fact, in both terms cr^J"* and c^J^, a fabric tensor can emerge 
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from the summation p!B|fT7|fT5] . The fabric tensor is useful to char- 
acterize the contact anisotropy of the assembly, which is known 
as a basic feature of granular assemblies. In dry granular materi- 
als, an induced anisotropy can develop when a deviatoric stress 
loading is applied. In partially saturated assembly, due to the pos- 
sibility of interactions without contact, the conclusion is not so 
trivial. As an illustration, if we restrict our analysis to spherical 
particles it can be inferred that: 

1 N N 

= 77 E E (8) 

y p=lq=l 

This relation points out that, contrary to the contact term where 
the fabric tensor ^ E^i nf ''^n^''^ is directly linked to the induced 
anisotropy (a contact is associated with a force), two causes can 
be invoked for the capillary term. First, the distribution of the 
liquid bonds can be anisotropic. Secondly, the geometry of the 
bonds being obviously dependent on the local geometry, it is 
possible that the distribution of both terms F^^p and /^'^ is also 
anisotropic. This is significant because the anisotropic attribute 
yields shear effects associated with pore fluid, which could mostly 
influence material behaviour. 

In order to enrich our discussion, we present numerical investiga- 
tions of these features in the following sections using both DEM 
and micromechanical simulations. 

3 MICROSTRUCTURAL INVESTIGATION ON THE CAPIL- 
LARY STRESS TENSOR USING DEM 

We present here a numerical analysis of the stress variables us- 
ing a micromechanical model based upon the Discrete Element 
Method initially introduced by Cundall and Strack ||2] . This tech- 
nique starts with basic constitutive laws between interacting par- 
ticles and can provide a macroscopic response of a particle assem- 
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bly due to loading changes at the boundaries. Each particle of the 
material is a rigid sphere identified by its own mass, m, radius, 
R and moment of inertia, Iq. For every time step of the com- 
putation, interaction forces between particles, and consequently 
resulting forces acting on each of them, are deduced from sphere 
positions through the interaction law. Newton's second law is 
then integrated through an explicit second order finite difference 
scheme to compute the new sphere positions. 

3.1 The discrete element model 

A 3D software called YADE (Yet Another Dynamic Engine), Koz- 
icki and Donze |[19] , has been enhanced in order to properly sim- 
ulate partially saturated granular material features. 

3.1.1 Inter-particle behaviour 

The contact interaction is described by an elastic-plastic relation 
between the force F and the relative displacement U of two in- 
teracting particles. A normal stiffness Kn is therefore defined to 
relate the normal force to the intergranular normal relative 
displacement Un '■ 

^ \ KnUn if Un<0 

Fn = (9) 

[ if Un>0 

and a tangential stiffness Kt allows us to deduce the shear force 
Ft induced by the incremental tangential relative displacement 
dUt'j this tangential behaviour obeying the Coulomb friction law: 

dFt = -KtdUt 

(10) 
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(a) D (b) 

Fig. 1. Illustration of a liquid bridge between two particles of unequal sizes: (a) 
global geometry, (b) details of the bridge. 

where fi is the Coulomb friction coefficient defined as /i = tan{(j))^ 
with (j) the intergranular friction angle. 



In the work presented here, and Kt are dependent functions 
of the interacting particle sizes and of a characteristic modulus 
of the material denoted as E: 

Kt = a Kr, 



This definition results in a constant ratio between E and the 
effective bulk modulus of the packing, whatever the size of the 
particles. 



For simplicity, we assume that the water inside the sample is 
solely composed of capillary water as defined in the pendular 
state, with a discontinuous liquid phase. 

Much attention has been given to these pendular liquid bridges. 
Their exact shape between spherical bodies is defined 



by the Laplace equation, relating the pressure difference Au = 
Ua — Uw across the liquid-gas interface to the mean curvature of 
the bridge and the surface tension of the liquid phase 7: 

Au = -f(- + -] (12) 
\ri r2/ 
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In the Cartesian coordinates of Fig|T](b) the two curvature radii 
ri and r2 (Fig|T](a)), are given by: 



- = , (13) 

ri y{x)Jl + y'^{x) 



and 

1 y"ix] 



T2 + 



(14) 



where y{x) defines the profile of the liquid-gas interface curve. 
The coincides with the axis of symmetry of the bridge, 

passing through the centers of the connected spheres (Fig|T](b)). 
According to the Laplace equation, the profile of the liquid bridge 
is thus related to the capillary pressure Am through the following 
differential equation: 

^(1 + y'\x)fl' + - y"{x) = (15) 



The corresponding liquid bridge volume V and intergranular dis- 
tance D can be obtained by considering the ^-coordinates {xi and 
X2) of the three-phases contact lines defining the solid-liquid-gas 
interface, as defined by Soulie et al., 



V = Ti y'^(^x)dx - \'nR\{l - acos{xi)Y{2 + acos{xi)) 

(16) 

— |7ri?2(l — acos(a:2))^(2 + acos{x2)) 



and 

D = i?2(l — acos{x2)) + X2 + Ri{l — acos{xi)) — xi (17) 

The capillary force due to the liquid bridge can be calculated 
at the profile apex yo according to the 'gorge method' [20] and 
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cap 




Fig. 2. Evolution of the capillary force Fcap with the intergranular distance D for 
a given suction value: a meniscus can form for D < Dcreation and breaks off for 

D > Dj-upture' 

consists of a contribution of the capillary pressure Au as well as 
the surface tension 7: 

Fcap = 27r?/o7 + ttz/qAw (18) 



The relation between the capillary pressure and the configuration 
of the capillary doublet is thus described by a system of non- 
linear coupled equations (15, 16, 17, 18) where the local geome- 
try [D) and water volume result of the solved system 
[ [23] . So, to account for capillarity in the code, an interpolation 
scheme on a set of discrete solutions of the Laplace equation has 
been developed in order to link directly the capillary pressure 
to the capillary force and water volume of the liquid bridge for a 
given grain-pair configuration. This results in a suction-controlled 
model where, at every time-step during the simulation, capillary 
forces and water volumes {Fcap, are computed based upon the 
microstructure geometry D and the imposed suction level Au. 

{Fcap,V) = ^{Au- D) (19) 



A schematic diagram of the implemented capillary law is shown 
in FigJ2]for a given value of the suction. In this paper, the choice 
was made to define the appearance of a bridge when particles 
come strictly in contact {Dcreation = 0), neglecting the possibility 
of adsorbed water effects. The capillary force is considered con- 
stant for the range of the elastic deformation {D < 0), assuming 
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local displacements to be very small compared to particle radii. 
Let us note that the formulation intrinsically defines the distance 
from which the meniscus breaks off as depending on the given 
capillary pressure and on the local geometry. This maximum dis- 
tance Drupture corrcspouds to the minimum D value from which 
the Laplace equation has no solution. 



3.1.2 Stress tensors 



Since this study covers the macroscopic and microscopic aspects 
of unsaturated granular media, stress tensors are calculated by 
both macro and micro-methods. The macro-method is the con- 
ventional way used in laboratory to measure stresses in experi- 
ments, that is to say: 

^^J = {T.F^)/S, (20) 



where Fj is the normal force acting on the boundary, and Sj the 
surface of the boundary oriented by the normal direction j. This 
stress is equivalent to net stress {aij — UaSij) used in unsaturated 
soils mechanics, with Ua used as the reference pressure because 
the pore air pressure is effectively zero in many practical prob- 
lems (as well as in this study). 



As seen in section 2.2, two other stress tensors can be considered 
through homogenization techniques: the intergranular stress ten- 
sor 0"'^°"* computed from intergranular forces, and the capillary 
stress tensor a^j^ computed from capillary forces. 

3.1.3 Sample description and testing programme 

The studied particle assembly is a 1 mm length cubic sample 
composed of 10000 spheres, with a grain size distribution ranging 
from 0.025 mm to 0.08 mm, as shown in FigJSl and a porosity 
of 0.385. The input parameters are listed in table 1, referring to 
equation (10). 
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Global Modulus 


Kt 


Friction angle 


E (MPa) 


a 


(j) (deg.) 


150 


0.5 


30 



Table 1 

DEM model parameters 




Grain sizes (mm) 



Fig. 3. Sample description. 

The sample was prepared by an isotropic-compaction technique, 
which can be described in two stages. 

(a) All particles are randomly positioned inside a cube made up 
of 6 rigid walls in such a manner that no overlap/contact force 
develops between any pair of particles. The interparticle fric- 
tion angle is set to a small value (the smaller the friction angle, 
the denser the assembly. A value of 0.2 degree has been cho- 
sen here) and particle radii are then homogeneously increased, 
whereas boundary walls stay fixed. The process is run until the 
confining pressure (15 kPa) is reached and equilibrium between 
the internal stress state and the external load is satisfied [21] • 

(b) The interparticle friction coefhcient is then changed to a value 
classically used in DEM simulation to reproduce an accept- 
able shear strength (30 degree) and boundary walls are servo- 
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controlled in displacement to keep the equilibrium state. 



Starting from the initially stable configuration, a given suction is 
then applied. In this case, capillary forces (Eq.(18)) are added to 
all existing contact forces {D creation = 0, see figure [2]). This pro- 
cess simulates the appearance of liquid bridges at body contacts 
as it would take place during a capillary condensation when the 
relative humidity of the surrounding air is increased, and coin- 
cides with a wetting of the material. The sample then reaches a 
new equilibium state from which different stress paths can be im- 
posed. Note that capillary forces are assumed to be zero between 
a sphere and a wall. 



Suction-controlled triaxial compression tests have been carried 
out on the generated specimen which was taken up from the ini- 
tial state of 15 kPa to upper confining pressures by isotropic com- 
paction through wall displacements. The capillary pressure value 
(10 /cPa) has been chosen so as liquid bridge volumes are small 
enough to avoid the possibility of interconnected liquid bridges, 
and hence to ensure the pendular regime in the medium. The 
water content can be simply computed as the sum of all the liq- 
uid bridge volumes. In this case, the associated initial degree of 
saturation of the sample is about 20 %. A constant compression 
rate is then applied in the axial direction, controlling the lateral 
walls in displacement to keep the confining pressure constant. 
To keep the quasi-static assumption, the loading rate was fixed 
sufficiently small so as the normalized mean resultant force on 
particles (which is at static equilibrium) does not exceed 1 % 
at each loading step. 
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Fig. 4. Axial stresses resulting from elastic forces (acont) and capillary forces {(Jcap)i 
as well as the external load on the wall (cr^aii) under a suction-controlled triaxial 
compression. 



3.2 Stress tensor analysis 



3.2.1 Macroscopic evidence for the capillary stress tensor 



The quasi-static assumption ensures the equihbrium of the in- 
ternal stress state and external load to be satisfied over every 
loading path. The internal stress has therefore to be the sum of 
the repulsive stress due to the elastic forces and the tensile 
one due to the capillary forces Fcap: as presented by Eq.(6) in 
subsection 2.2. 



As shown in figure [4] for a triaxial loading path under a given 
confining pressure of 15kPa, this additivity is perfectly verified. 
The quasi-static assumption is numerically confirmed and the 
existence of a capillary stress tensor is therefore proved for the 
case of wet granular materials. 




cap 
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Fig. 5. Evolution of the principal capillary stress tensor components during a suc- 
tion-controlled triaxial compression test {Pq = IbkPa, u = WkPa). 

3.2.2 Capillary stress tensor analysis 

This suction induced stress leads to some questions about its 
structure and more generally on water distribution inside the 
material. In fact, classical considerations of unsaturated mate- 
rials often assimilate suction effects to an equivalent pressure 
which consequently acts in the medium independently from its 
anisotropy (the "hydraulic" component of the effective stress is 
generally considered as an isotropic tensor). 



Computing the principal components of this capillary stress ten- 
sor along the deviatoric loading path of a triaxial test compro- 
mises this assumption as shown in Figj5l This is all the more 
remarkable in that the model ensures a uniform distribution of 
the capillary pressure inside the medium. It is clear that for the 
initial state corresponding to an isotropic configuration of the 
assembly, the capillary stress tensor is almost spherical with an 
initial mean value of about 3.6 kPa for both axial and lateral 
components. Nevertheless, the anisotropy rapidly evolves with 
the one induced by loading, providing a difference between the 

cap cap 

principal components. If we define a = 2 l^lap^lup^ as a represen- 
tative index of the tensor sphericity: a = for the initial isotropic 
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epsl {%} epsl (%) 

Fig. 6. Evolutions of the porosity and average number of contacts by particles K 
under a suction-controlled triaxial compression {Pq = IbkPa, u = 10 kPa). 

state, and then slightly evolves to a quasi-constant value of 0.12 
from the 7 % deformation level until the final 15 % computed 
state. 

The causes of this evolution can be analysed by examining the 
volumetric deformation of the sample (FigJSl(a)) and the associ- 
ated packing rearrangement through the average number of con- 
tacts by particle, K, (Figi6](b)). As these global considerations 
have to be completed by a micromechanical analysis in order to 
gain a clear insight into the microstructural origins of the phe- 
nomenon, FigiG] will be commented in the following section. 

3.2.3 Micromechanical investigation 

Here we develop a micromechanical analysis of the stress vari- 
ables by considering both contact and liquid bridge distributions 
through the whole assembly. As the local interactions (dry con- 
tacts or menisci interactions) involve normal directions, a database 
can be defined in terms of orientations for all the grains of the 
sample. It is proposed to examine this normal direction network. 

The search has been done considering the given direction angle 
as presented in FigJH with 6 corresponding to the angle of 
the unit normal vector (n) from the axis of axisymmetry of the 
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Fig. 7. Description of the method used for interaction orientation distribution. 

sample (Y). As seen in section 2.2, a^J^ can be written as: 

1 N N 

= ^ E E J'S'''''"f"'nr' (22) 

V p=lq=l 

which points out the possible induced anisotropy of the capil- 
lary tensor by means of both liquid bridge and force intensity 
distributions. If we now introduce Pmeniscus{n) as the menisci ori- 
entation distribution inside the sample (in fact, the number of 
menisci along the direction n) defined by Iv Pmeniscus{n)dV = 1, 
Eq.(22) becomes: 

I r 

(T^j^ = — J < Fcap.l >n Pmeniscus{n)n ® fidV (23) 
V 

— > — # 

where < F^ap-l >n is the mean value of F^ap-l along the direc- 
tion n. It is therefore possible to compute separately the geomet- 
ric distribution {Pmeniscus{n)) and the static distribution of the 
(< Fcap-l >n) quantity which involves the mean force intensity, 
for every direction characterized by n. 

P{n) (FigiHj) and < F.l >^ (fig-E]) distributions with 9 are plotted 
for several deformation levels for both contacts and menisci con- 
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Fig. 8. Contacts and menisci orientation distribution -P(n) for different deformation 
levels 

tributions. < Fcojit-l >n and < Fcap-l >n are simply normalized 
by their mean value to be qualitatively compared. The search 
has been done for different deformation levels on the sample con- 
fined under 15 kPa and subjected to a constant capillary pressure 
of 10 kPa. Different snapshots have been taken, starting from 
the assumed isotropic initial state, until a 15 % axial strain level 
where the deformation regime appears almost permanent. First, 
as menisci are added at contacts to simulate capillary condensa- 
tion, distributions of both contact and capillary terms are iden- 
tical in the initial state. The structural isotropy of the sample 
clearly appears with a uniform distribution of P{n) and < /./ >^ 
for contacts and menisci along all the directions (Figsj8](a) and 
[9](a)), confirming here the accuracy of the generation process. 
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The evolutions of both contact and hquid bridge distributions 
during the loading directly results from the deformations of the 
assembly. As the sample reacts just like a dense granular material 
(FigJ5](a)), the initial contractancy gives rise to a brief increase in 
the coordination number (FigJSl(b)), leading to the development 
of new liquid bridges. The corresponding growth of both axial 
and lateral capillary stress components can be viewed on FigJS 
As a consequence of the persistance of menisci for low interparti- 
cle distances, this combined augmentation persists until ei = 3 % 
even if K strongly drops before rising up again. The small dif- 
ference between the axial and the lateral components is simply 
due to the deviatoric loading which produces more contacts, and 
consequently more menisci, in the active loading direction (Y) 
than in the passive stress-controlled one. 

After ei = 3%, the lateral capillary tensor component clearly 
starts to reduce. As pointed out in FigJH](b), this results from 
the lateral spreading of the particles produced by the dilatancy 
of the assembly. Even though menisci distribution Pmeniscus{n) 
seems not to follow the induced fabric anisotropy because of the 
remaining of the liquid bridges, lateral capillary forces tend to 
diminish due to increasing interparticular distances (Fig|2]). On 
the other hand, the axial component of the capillary stress ten- 
sor regularly rises until the initiation of a permanent regime in 
the deformation process (near ei = 7%) where K continues on 
going quasi-constant. From this state, as the number of contacts 
stabilizes while dilatancy persists, < F^ap-l >n cannot endure any 
further increase and the axial component of the capillary tensor 
starts to diminish because of a higher spreading of the grains. 

Describing the fabric induced anisotropy, Pcont{n) distribution 
from 3 to 15%, seems to be quite constant. However, it can be 
noted that < Fcont-l >n contribution appears to reach rapidly a 
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(c) epsl = 7% (d) epsl = 15% 

Fig. 9. Contacts and capillary forces distribution orientations (< F.l >fi) for differ- 
ent deformation levels 

maximum anisotropic state before slightly reducing to the final 
one. This maximum anisotropic strength state fairly corresponds 
to the peak shear strength of the sample (ei = 3 %) where verti- 
cal contact force chains are subjected to maximum loading before 
breaking off. Concerning the final distributions, they are well rep- 
resentative of the so-called critical state, where a stabilization of 
the stresses occurs with no new significative change in the evolu- 
tion of the anisotropy. In a noteworthy way, menisci seem not to 
follow the same evolution. Effectively, as liquid bridges can ex- 
ist in a certain range of increasing intergranular distances, their 
distribution in the media is not driven in the same way by the fab- 
ric induced anisotropy and tends to stay close to the initial state, 
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particularly for the range of small deformations. Nevertheless, for 
large deformations, due to the sample dilatancy, a small induced 
anisotropy arises from the disappearance of liquid bridges in the 
lateral directions. It is evident that the history of the material is 
fundamentally essential when dealing with water distribution. 



To sum up, the analysis reveals a slight induced anisotropy of 
the capillary stress tensor as a function of the medium fabric. 
The pore fluid in unsaturated soil has its own fabric that may 
be readily altered with changes in the granular fabric and is also 
strongly dependent on the water distribution inside the media. 
The global approximation which characterizes water effects in 
unsaturated materials by an equivalent pore pressure is, there- 
fore, unable in essence to point out this intrinsically anisotropic 
microstructural force contribution. However, depending on the 
hydric history, the evolving anisotropy of the pore water distri- 
bution can validate the assumption by counterbalancing the in- 
duced fabric anisotropy in the material. 



Since DEM analyses result from direct simulations of a granular 
assembly, the purpose of the next section is to compare the re- 
sults with those of a microstructural model where the behaviour 
of the material is obtained through a micromechanically based 
constitutive relation. 



4 COMPARISON WITH A MICRO-MECHANICAL MODEL 



In this section, we first present the microstructural model used 
to compare with DEM simulations. This is a stress-strain model 
( |25p6l fT]) proposed by Chang and Hicher which considers 
inter-particle forces and displacements. Its capability has been 
recently extended to unsaturated states by incorporating the in- 
fluence of capillary forces at the micro level jlO]. By comparing 
the predicted triaxial loading results obtained by the two ap- 
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proaches for the granular assembly of section 3, we confirm the 
stress conceptions introduced previously, focusing on the capil- 
lary stress tensor as defined in section 2.2. 

4.1 Stress-Strain Model 

In this model, we envision a granular material as a collection 
of particles. The deformation of a representative volume of this 
material is generated by the mobilisation of contact particles in 
all orientations. Thus, the stress-strain relationship can be de- 
rived as an average of the mobilisation behaviour of local contact 
planes in all orientations. The forces and movements at the con- 
tact planes of all orientations are suitably superimposed to obtain 
the macroscopic stress strain tensors using the static homogeniza- 
tion presented in Section 2.2. 

4-1.1 Inter-particle behaviour 

- Contact forces and capillary forces 

For dry samples, contact forces are directly determined from the 
external stresses a applied on the granular assembly (Eq.(2)). In 
the case of wetted samples, different stages of saturation can be 
identified. The fully saturated regime corresponds to a two-phase 
material with water filling completely the voids between grains. 
The water pressure Uy^ can either be positive or negative (suction) 
but in both cases the effective stress concept [1] can be applied 
and the contact forces determined by considering the effective 
stresses cr' as the external stresses ([ I28.29J ): 

a = a- u^I (24) 

As seen previously, in the case of partially saturated samples in 
the pendular regime, the liquid phase is distributed in menisci lo- 
cated between close grains. As a consequence, capillary forces are 
applied on the grains and are added to the contact forces defined 
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above. The attractive capillary force between two grains con- 
nected by a water bridge is a decreasing function of the distance 
between the grains until the bridge fails (FigJ2])- This function de- 
pends on the volume of liquid found between the grains. Different 
mathematical expressions have been proposed for these capillary 
forces, Fcap- Eq.(18) presents the expression used in DEM. F^ap 
depends on the capillary pressure defined as the pressure jump 
across the liquid-air interface, on the liquid-air interface surface 
tension, as well as on the geometry of the menisci governed by the 
solid-liquid contact angle and the filling angle. One can see that 
Fcap depends on the geometry of the liquid bridge which is func- 
tion of the amount of the pore water and of the distance between 
two neighboring grains. The use of Eq.(18) for determining the 
amplitude of the capillary forces is not a straightforward one and 
therefore, in the micro-mechanical model, a simplified approach 
is to consider an empirical relation between F^ax and the degree 
of saturation Sr: 

F"^ = F e'^^^^ (25) 

^ cap ^ max^ K'^^J 



where F^^^ is the capillary force between two neighbouring grains, 
not necessarily in contact, Fmax the value of F^ap for two grains in 
contact and R the mean grain radius. D represents the distance 
between two grains and is equal to / — 2R, I being the branch 
length given as a distribution function of the grain size and the 
void ratio and c is a material parameter, dependent on the grain 
morphology and on the water content: 



Fmax — Fqj^ for < Sr < Sq 
Fmax = -^O 5'°(i-5'o) ^r Sq < 5,. < 



(26) 



where Fq and Sq are material parameters. Fq depends on the 
grain size distribution, 5*0 represents the degree of saturation at 
which any further drying of the specimen will cause substantial 
breaking of the menisci in the pendular domain. 5*0 depends on 
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the nature of the granular material. 



Since the menisci are not necessarily all formed in the funicular 
regime, Eq.(25) may not be applicable for high degrees of satura- 
tion. However, in this first approach, we decided to extend it to 
the whole range of saturation, considering that the amplitudes of 
capillary forces were small for degrees of saturation higher that 
80% and could thus be approached with sufficient accuracy by 
using the same equation. 

- Elastic relationship 

The contact stiffness of a contact plane includes normal stiffness, 
/c" , and shear stiffness, k^. The elastic stiffness tensor is defined 
by 

F^^ = A;f/(5f (27) 

which can be related to the contact normal and shear stiffness by 

kff = k»] + k^s^s'^ + tft]) (28) 

The value of the stiffness for two elastic spheres can be estimated 
from Hertz-Mindlin's fomulation. For sand grains, a revised form 
was adopted pD], given by 

A:n = U(^)" k = ho{^,r (29) 

where Gg is the elastic modulus for the grains, F„ is the contact 
force in normal direction, / is the branch length between the two 
particles, kno, ho and n are material constants. 

- Plastic relationship 

Plastic sliding often occurs along the tangential direction of the 
contact plane with an upward or downward movement, thus shear 
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dilation/contraction takes place. The dilatancy effect can be de- 
scribed by 

^ = ^ - *»0„ (30) 



wfiere 0o is a material constant which, in most cases, can be con- 
sidered equal to the internal friction angle 0^. This equation can 
be derived by equating the dissipation work due to plastic move- 
ments and friction in the same orientation. Note that the shear 
force T and the rate of plastic sliding dA^ can be defined as 
T = /fJTW and d/V = ^{d5^Y + {dd^f. The yield func- 
tion is assumed to be of the Mohr-Coulomb type, 

^) = T- F^AP) = (31) 



where k{Ap) is an isotropic hardening/softening parameter de- 
fined as: 

^ ^ kpotan{(l)p)AP ^^^^ 
\Fn\tan{(J)p) + kpoAP 



The hardening function is defined by a hyperbolic curve in k — A^ 
plane, which involves two material constants: (f)p and KpQ. On the 
yield surface, under a loading condition, the shear plastic flow 
is determined by a normality rule applied to the yield function. 
However, the plastic flow in the direction normal to the contact 
plane is governed by the stress-dilatancy equation in Eq.(32). 
Thus, the flow rule is non-associated. 

- Interlocking influence 

The internal friction angle 0^ is a constant for the material. How- 
ever, the peak friction angle, 0^, on a contact plane is dependent 
on the degree of interlocking by neighboring particles, which can 
be related to the state of the packing void ratio e by: 

tan{(f)p) = {-rtan{(j)^) (33) 
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where m is a material constant [21] • The state of packing is itself 
related to the void ratio at critical state Cc- The critical void 
ratio Cc is a function of the mean stress. The relationship has 
traditionally been written as: 

P' 

Cc = r - \log{p') or ec = e^g/ - A/o^( ) (34) 

Pref 

where T and A are two material constants and p' is the mean 
stress of the packing, and (e^e/, Pref) is a reference point on the 
critical state line. 

For dense packing, the peak frictional angle 0^ is greater than (f)^. 
When the packing structure dilates, the degree of interlocking and 
the peak frictional angle are reduced, which results in a strain- 
softening phenomenon. 

- Elasto-plastic relationship 

With the elements discussed above, the final incremental stress- 
strain relations of the material can be derived that includes both 
elastic and plastic behaviour, given by = k'^fS'j. Detailed 
expression of the elasto-plastic stiffness tensor is given in [27] . 



4.1.2 Stress-strain relationship 
- Macro micro relationship 

The stress-strain relationship for an assembly can be determined 
from integrating the behaviour of inter-particle contacts in all 
orientations. During the integration process, a relationship is re- 
quired to link the macro and micro variables. Using the static 
hypotheses proposed by Liao et. al [32], we obtain the relation 
between the macro strain and inter-particle displacement (finite 
strain condition not being considered here): 

N 

= Ar,^ E S^lk (35) 

a=l 
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where is the relative displacement between two contact parti- 
cles and the branch vector Ik is the vector joining the centers of 
two contacting particles. 

Using both the principle of energy balance and Eq.(36), the mean 
force on the contact plane of each orientation is 

= ^^jAjk'lkV (36) 

The stress increment aij induced by the loading can then be ob- 
tained through the contact forces and branch vectors for contacts 
in all orientations pTjfTB] . Since 1^ represents the mean branch 
vector for the a^'^ orientation including both contact and non- 
contact particles, the value of in Eq.(37) represents the mean 
of contact forces in the a^^ orientation. 

= ^ E F^l] (37) 

When the defined contact force is applied in Eq.(37), Eq.(38) is 
unconditionally satisfied. 



Using the definition of Eq.(38), the stress induced by capillary 
forces can be computed and is termed as capillary stress, given 
by 

cap 1 ^ a a 
^ij ~ 77 E F^ap,i^'j (38) 
^ a=l 

As mentionned in section 2.2, it is noted that this term is not anal- 
ogous to the usual concept of capillary pressure or suction which 
represents the negative pore water pressure inside the unsatu- 
rated material. In agreement with the results obtained through 
DEM simulations, the capillary stress depends on the geometry 
of the pores and is a tensor rather than a scalar. Only for an 
isotropic distribution of the branch lengths this water associ- 
ated stress can be reduced to an isotropic tensor. This is the case 
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for an initially isotropic structure during isotropic loading, but 
during deviatoric loading, an induced anisotropy is created and 
the capillary tensor is no longer isotropic. 



At the equilibrium state, the effective intergranular forces will 
therefore be the difference between the repulsive forces due to 
the external stresses (Eq.(38)) and the attractive capillary forces 
(Eq.(39)). In a similar way to Eq.(7), we can thus define a gen- 
eralized intergranular stress tensor a*, defined by: 

a' = a- a'^P (39) 



Assuming that a* can stand as an appropriate definition of the 
effective stress, this equation could represent a generalisation of 
Eq.(l) proposed by Bishop, in which the capillary stress is re- 
duced to an isotropic tensor. Dangla et al. ||33] demonstrated the 
validity of the effective stress approach in elasticity by means 
of an energy approach. They obtained an expression similar to 
Eq.(40) but with an additional term corresponding to the work 
of the interfaces. As pointed out before, capillary forces in our 
models, and consequently the capillary stresses, depend on the 
negative pore water pressure, or suction, and on the water-air 
interface surface tension. A similar approach can be found in the 
work of Fleureau and al. [M] who obtained an explicit expres- 
sion of the capillary stress as a function of the suction for regular 
arrangements of spherical grains by neglecting the surface ten- 
sion. The definition of the capillary stress tensor in Eq.(39) can 
therefore be considered as an extension of the results obtained 
from these previous studies to cases of non isotropic granular 
assemblies. 

4- 1-3 Summary of parameters 

One can summarize the material parameters as: 
- Normalized contact number per unit volume: 4?-. 
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- mean particle size, 2R. 

- Inter-particle elastic constants: /c^o? and n. 

- Inter-particle friction angle: (j)^ and m. 

- Inter-particle hardening rule: kpQ and 0o- 

- Critical state for packing: A and F or e^e/ and Pref- 

- Capillary force equation: /q, Sq and c. 

Other than critical state parameters, all parameters are inter- 
particles. Standard values for kpQ and 0o are the following: kpQ = 
/c„o and 0o = 0y and a typical ratio = 0.4 can generally 
be assumed jTD]. Therefore, for dry or saturated samples, only 
six parameters have to come from experimental results and these 
can all be determined from the stress-strain curves obtained from 
drained or undrained compression triaxial tests. For unsaturated 
states, three more parameters need to be determined, using spe- 
cific triaxial tests on partially saturated samples. 

4-2 Numerical simulations 

Several simulations of triaxial loading have been performed in or- 
der to compare DEM and Micromechanical Model results. These 
are based on computations on a Representative Volume Element 
of about 10 000 spherical elements. 

4-2.1 Dry samples 

The model parameters were determined from the following pro- 
cedure. The granular assembly defined in the discrete element 
model is made of spherical particles with the grain size distribu- 
tion presented in FigJ21 The particle size 2R was selected equal 
to = 0.045 mm. The elastic parameters could not be directly 
derived from the inter-particle behaviour used in the DEM sim- 
ulations, which consider a linear contact stiffness. From previ- 
ous studies on glass beads assemblies pSlfTD] , typical values were 
adopted in this study. The plastic parameters were determined 
from the numerical results obtained by DEM. Table II summa- 
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rizes the set of parameters used for modelling the dry sample 
behaviour. 



kno (N/m) 


#2. n 


(pfj, (deg.) 


A m 


300 


0.5 30 


20 


0.05 0.5 



Table 2 

Microstructural Model parameters for the glass beads assembly 

FigHO] presents the numerical simulations of three triaxial tests 
performed at three different confining pressures of 15, 30 and 
60 kPa for the simulated dry assembly having an initial void ra- 
tio equal to 0.38. One can see that the results obtained with the 
Microstructural Model compared well to DEM ones. 

4-2.2 Wet samples 

DEM simulations were performed on unsaturated assemblies with 
an initial saturation degree of about 20 %, corresponding to a suc- 
tion value equal to 10 kPa. In order to determine the correspond- 
ing capillary forces, Eq.(25) includes two material parameters c 
and d which control the evolution of the water induced forces with 
the distance between two particles. According to experimental re- 
sults presented in different studies |2Tll23j . we decided to take a 
value of c = 4 and of d = 0.05 mm. These values give a standard 
evolution for the capillary forces, function of the distance between 
particles, as well as an initial isotropic distribution of these forces 
if the material structure is isotropic. The evolution of the capil- 
lary forces with the degree of saturation requires two more param- 
eters: 5*0 and /q. According to previous studies |ID], we selected 
a value of 5*0 = 1 % and determined the value of /o = 0.045 N. 
In order to obtain an initial value of the capillary stress in ac- 
cordance with the one computed by DEM. We then performed 
numerical simulations of triaxial tests on wet samples for dif- 
ferent confining pressures. Contrary to DEM suction-controlled 
simulations, the Microstructural Model tests are water content 
controlled. In order to compare those two approaches, samples 
were initially wetted at a common degree of saturation of about 
20 %. One should notice that, even if the test conditions are not 
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Fig. 10. DEM and MicroMechanical simulations of triaxial compression tests on a 
similar dry granular assembly. 

strictly identical, the changes in the degree of saturation dur- 
ing loading obtained for both tests are sufficiently similar for us 
(FiglH]) to compare the results obtained by the two approaches. 



As presented in FigHSl the two models give quite similar results. 
One can see that a material strength increase is obtained for un- 
saturated samples compared to dry ones at the same confining 
pressure. The volume changes during triaxial loading create a 
small change in the degree of saturation (FiglH]). As a conse- 
quence, the capillary forces evolve, according to Eq.(5). FigHSl 
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Fig. 11. Evolution of the saturation degree during DEM and Micromechanical sim- 
ulated triaxial tests. 

shows the evolution of the principal components of the capillary 
stress tensor during constant water content triaxial tests. The 
initial state corresponds to an isotropic capillary stress tensor 
with a mean stress equal to 3.6 kPa as obtained in DEM. During 
loading, a structural anisotropy is created due to the evolution 
of the fabric tensor defined in Eq.(38). Therefore, the principal 
components of the capillary stress tensor evolve differently. In 
the studied cases, this difference remains small and corresponds 
at the end of the test to a relative difference less than 10 %. This 
small difference can be explained by two causes. The first one is 
the small amount of induced anisotropy obtained by the evolution 
of the fabric tensor. This evolution is due to the change in the 
branch length If for each a direction which, in this version of the 
model, is only due to elastic deformations of the grains in contact. 
Since all our numerical testing were performed at small confin- 
ing stresses, the amount of elastic deformation remained quite 
small. The second reason is linked to the fact that the capillary 
bridges can exist for non-touching neighboring grains. This has 
been taken into account in calculating the mean capillary force 
and also in determining the capillary stress tensor (Eqs.(27) and 
(41)). This result is in agreement with the distribution of the 
contacts and menisci distribution computed by DEM (FigsJHland 

El). 
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Fig. 12. DEM and Micromechanical simulations of triaxial compression tests on a 
similar wet granular assembly. 



Regarding the constitutive behaviour at contacts between solid 
particles, the results provided by the two methods are in rather 
good agreement concerning the influence of capillary forces at a 
macroscopic level. The increase in the shear strength classically 
observed for partially saturated materials is clearly encountered 
starting from microscopical considerations, and the slight induced 
anisotropy of the capillary stress tensor is observed, confirming 
that suction effects in unsaturated materials cannot be precisely 
accounted for by an equivalent pore pressure assumption. 
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Fig. 13. Evolution of the principal capillary stress tensor components during a tri- 
axial compression test on a wet granular assembly. 

5 CONCLUSION 



Starting from local capillary forces, a stress variable, denoted 
as the capillary stress tensor and intrinsically associated to wa- 
ter effects, has been defined through homogenisation techniques. 
Triaxial compression test simulations from two fundamentally 
different micromechanical models were performed on a granu- 
lar assembly under several confining pressures for dry and par- 
tially saturated conditions. Both models reproduce in quite good 
agreement the main features of unsaturated granular materials, 
in particular the increase of the shear strength due to capillary 
effects. 

The results also suggest that, in partially saturated materials 
within the pendular regime, the effects of pore fluid are ade- 
quately represented by a discrete distribution of forces rather 
than by an averaged pressure in the fluid phase. Effectively, as a 
representative quantity of the pore fluid distribution inside unsat- 
urated materials, this suction associated stress tensor reveals that 
pore fluid has its own fabric which is inherently anisotropic and 
strongly dependent on the combined loading and hydric history 
of the material. Even if the induced anisotropy of the capillary 
stress tensor appears slight in this study, it is obvious that this 
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tensorial nature of water in unsaturated material implies suc- 
tion to produce shear effects on the solid phase. This suction in- 
duced shear effect consequently makes it difficult to associate an 
isotropic quantity to water as expressed in the Bishop's effective 
stress. Pore pressure is no longer an isotropic stress in unsatu- 
rated soil, and therefore, cannot be considered as an equivalent 
continuum medium. The analysis finally confirms that suction is 
a pore-scale concept, and that stress definitions for unsaturated 
soils should also include microscopic interparticle stresses as the 
ones resulting from capillary forces. 

The multi-scale approach presented here appears to be a perti- 
nent complementary tool for the study of unsaturated soil me- 
chanics. More precisely, discrete methods should convey new in- 
sights into the discussion about the controversial concept of gen- 
eralized effective stress by relating basic physical aspects to clas- 
sical phenomenological views. 
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